home *** CD-ROM | disk | FTP | other *** search
/ IRIX Base Documentation 2002 November / SGI IRIX Base Documentation 2002 November.iso / usr / share / catman / p_man / cat3 / SCSL / stgsna.z / stgsna
Encoding:
Text File  |  2002-10-03  |  14.6 KB  |  397 lines

  1.  
  2.  
  3.  
  4. SSSSTTTTGGGGSSSSNNNNAAAA((((3333SSSS))))                                                          SSSSTTTTGGGGSSSSNNNNAAAA((((3333SSSS))))
  5.  
  6.  
  7.  
  8. NNNNAAAAMMMMEEEE
  9.      STGSNA - estimate reciprocal condition numbers for specified eigenvalues
  10.      and/or eigenvectors of a matrix pair (A, B) in generalized real Schur
  11.      canonical form (or of any matrix pair (Q*A*Z', Q*B*Z') with orthogonal
  12.      matrices Q and Z, where Z' denotes the transpose of Z
  13.  
  14. SSSSYYYYNNNNOOOOPPPPSSSSIIIISSSS
  15.      SUBROUTINE STGSNA( JOB, HOWMNY, SELECT, N, A, LDA, B, LDB, VL, LDVL, VR,
  16.                         LDVR, S, DIF, MM, M, WORK, LWORK, IWORK, INFO )
  17.  
  18.          CHARACTER      HOWMNY, JOB
  19.  
  20.          INTEGER        INFO, LDA, LDB, LDVL, LDVR, LWORK, M, MM, N
  21.  
  22.          LOGICAL        SELECT( * )
  23.  
  24.          INTEGER        IWORK( * )
  25.  
  26.          REAL           A( LDA, * ), B( LDB, * ), DIF( * ), S( * ), VL( LDVL,
  27.                         * ), VR( LDVR, * ), WORK( * )
  28.  
  29. IIIIMMMMPPPPLLLLEEEEMMMMEEEENNNNTTTTAAAATTTTIIIIOOOONNNN
  30.      These routines are part of the SCSL Scientific Library and can be loaded
  31.      using either the -lscs or the -lscs_mp option.  The -lscs_mp option
  32.      directs the linker to use the multi-processor version of the library.
  33.  
  34.      When linking to SCSL with -lscs or -lscs_mp, the default integer size is
  35.      4 bytes (32 bits). Another version of SCSL is available in which integers
  36.      are 8 bytes (64 bits).  This version allows the user access to larger
  37.      memory sizes and helps when porting legacy Cray codes.  It can be loaded
  38.      by using the -lscs_i8 option or the -lscs_i8_mp option. A program may use
  39.      only one of the two versions; 4-byte integer and 8-byte integer library
  40.      calls cannot be mixed.
  41.  
  42. PPPPUUUURRRRPPPPOOOOSSSSEEEE
  43.      STGSNA estimates reciprocal condition numbers for specified eigenvalues
  44.      and/or eigenvectors of a matrix pair (A, B) in generalized real Schur
  45.      canonical form (or of any matrix pair (Q*A*Z', Q*B*Z') with orthogonal
  46.      matrices Q and Z, where Z' denotes the transpose of Z. (A, B) must be in
  47.      generalized real Schur form (as returned by SGGES), i.e. A is block upper
  48.      triangular with 1-by-1 and 2-by-2 diagonal blocks. B is upper triangular.
  49.  
  50.  
  51.  
  52. AAAARRRRGGGGUUUUMMMMEEEENNNNTTTTSSSS
  53.      JOB     (input) CHARACTER*1
  54.              Specifies whether condition numbers are required for eigenvalues
  55.              (S) or eigenvectors (DIF):
  56.              = 'E': for eigenvalues only (S);
  57.              = 'V': for eigenvectors only (DIF);
  58.              = 'B': for both eigenvalues and eigenvectors (S and DIF).
  59.  
  60.  
  61.  
  62.  
  63.                                                                         PPPPaaaaggggeeee 1111
  64.  
  65.  
  66.  
  67.  
  68.  
  69.  
  70. SSSSTTTTGGGGSSSSNNNNAAAA((((3333SSSS))))                                                          SSSSTTTTGGGGSSSSNNNNAAAA((((3333SSSS))))
  71.  
  72.  
  73.  
  74.      HOWMNY  (input) CHARACTER*1
  75.              = 'A': compute condition numbers for all eigenpairs;
  76.              = 'S': compute condition numbers for selected eigenpairs
  77.              specified by the array SELECT.
  78.  
  79.      SELECT  (input) LOGICAL array, dimension (N)
  80.              If HOWMNY = 'S', SELECT specifies the eigenpairs for which
  81.              condition numbers are required. To select condition numbers for
  82.              the eigenpair corresponding to a real eigenvalue w(j), SELECT(j)
  83.              must be set to .TRUE.. To select condition numbers corresponding
  84.              to a complex conjugate pair of eigenvalues w(j) and w(j+1),
  85.              either SELECT(j) or SELECT(j+1) or both, must be set to .TRUE..
  86.              If HOWMNY = 'A', SELECT is not referenced.
  87.  
  88.      N       (input) INTEGER
  89.              The order of the square matrix pair (A, B). N >= 0.
  90.  
  91.      A       (input) REAL array, dimension (LDA,N)
  92.              The upper quasi-triangular matrix A in the pair (A,B).
  93.  
  94.      LDA     (input) INTEGER
  95.              The leading dimension of the array A. LDA >= max(1,N).
  96.  
  97.      B       (input) REAL array, dimension (LDB,N)
  98.              The upper triangular matrix B in the pair (A,B).
  99.  
  100.      LDB     (input) INTEGER
  101.              The leading dimension of the array B. LDB >= max(1,N).
  102.  
  103.      VL      (input) REAL array, dimension (LDVL,M)
  104.              If JOB = 'E' or 'B', VL must contain left eigenvectors of (A, B),
  105.              corresponding to the eigenpairs specified by HOWMNY and SELECT.
  106.              The eigenvectors must be stored in consecutive columns of VL, as
  107.              returned by STGEVC.  If JOB = 'V', VL is not referenced.
  108.  
  109.      LDVL    (input) INTEGER
  110.              The leading dimension of the array VL. LDVL >= 1.  If JOB = 'E'
  111.              or 'B', LDVL >= N.
  112.  
  113.      VR      (input) REAL array, dimension (LDVR,M)
  114.              If JOB = 'E' or 'B', VR must contain right eigenvectors of (A,
  115.              B), corresponding to the eigenpairs specified by HOWMNY and
  116.              SELECT. The eigenvectors must be stored in consecutive columns ov
  117.              VR, as returned by STGEVC.  If JOB = 'V', VR is not referenced.
  118.  
  119.      LDVR    (input) INTEGER
  120.              The leading dimension of the array VR. LDVR >= 1.  If JOB = 'E'
  121.              or 'B', LDVR >= N.
  122.  
  123.      S       (output) REAL array, dimension (MM)
  124.              If JOB = 'E' or 'B', the reciprocal condition numbers of the
  125.              selected eigenvalues, stored in consecutive elements of the
  126.  
  127.  
  128.  
  129.                                                                         PPPPaaaaggggeeee 2222
  130.  
  131.  
  132.  
  133.  
  134.  
  135.  
  136. SSSSTTTTGGGGSSSSNNNNAAAA((((3333SSSS))))                                                          SSSSTTTTGGGGSSSSNNNNAAAA((((3333SSSS))))
  137.  
  138.  
  139.  
  140.              array. For a complex conjugate pair of eigenvalues two
  141.              consecutive elements of S are set to the same value. Thus S(j),
  142.              DIF(j), and the j-th columns of VL and VR all correspond to the
  143.              same eigenpair (but not in general the j-th eigenpair, unless all
  144.              eigenpairs are selected).  If JOB = 'V', S is not referenced.
  145.  
  146.      DIF     (output) REAL array, dimension (MM)
  147.              If JOB = 'V' or 'B', the estimated reciprocal condition numbers
  148.              of the selected eigenvectors, stored in consecutive elements of
  149.              the array. For a complex eigenvector two consecutive elements of
  150.              DIF are set to the same value. If the eigenvalues cannot be
  151.              reordered to compute DIF(j), DIF(j) is set to 0; this can only
  152.              occur when the true value would be very small anyway.  If JOB =
  153.              'E', DIF is not referenced.
  154.  
  155.      MM      (input) INTEGER
  156.              The number of elements in the arrays S and DIF. MM >= M.
  157.  
  158.      M       (output) INTEGER
  159.              The number of elements of the arrays S and DIF used to store the
  160.              specified condition numbers; for each selected real eigenvalue
  161.              one element is used, and for each selected complex conjugate pair
  162.              of eigenvalues, two elements are used.  If HOWMNY = 'A', M is set
  163.              to N.
  164.  
  165.      WORK    (workspace/output) REAL array, dimension (LWORK)
  166.              If JOB = 'E', WORK is not referenced.  Otherwise, on exit, if
  167.              INFO = 0, WORK(1) returns the optimal LWORK.
  168.  
  169.      LWORK   (input) INTEGER
  170.              The dimension of the array WORK. LWORK >= N.  If JOB = 'V' or 'B'
  171.              LWORK >= 2*N*(N+2)+16.
  172.  
  173.              If LWORK = -1, then a workspace query is assumed; the routine
  174.              only calculates the optimal size of the WORK array, returns this
  175.              value as the first entry of the WORK array, and no error message
  176.              related to LWORK is issued by XERBLA.
  177.  
  178.      IWORK   (workspace) INTEGER array, dimension (N + 6)
  179.              If JOB = 'E', IWORK is not referenced.
  180.  
  181.      INFO    (output) INTEGER
  182.              =0: Successful exit
  183.              <0: If INFO = -i, the i-th argument had an illegal value
  184.  
  185. FFFFUUUURRRRTTTTHHHHEEEERRRR DDDDEEEETTTTAAAAIIIILLLLSSSS
  186.      The reciprocal of the condition number of a generalized eigenvalue w =
  187.      (a, b) is defined as
  188.  
  189.           S(w) = (|u'Av|**2 + |u'Bv|**2)**(1/2) / (norm(u)*norm(v))
  190.  
  191.      where u and v are the left and right eigenvectors of (A, B) corresponding
  192.  
  193.  
  194.  
  195.                                                                         PPPPaaaaggggeeee 3333
  196.  
  197.  
  198.  
  199.  
  200.  
  201.  
  202. SSSSTTTTGGGGSSSSNNNNAAAA((((3333SSSS))))                                                          SSSSTTTTGGGGSSSSNNNNAAAA((((3333SSSS))))
  203.  
  204.  
  205.  
  206.      to w; |z| denotes the absolute value of the complex number, and norm(u)
  207.      denotes the 2-norm of the vector u.
  208.      The pair (a, b) corresponds to an eigenvalue w = a/b (= u'Av/u'Bv) of the
  209.      matrix pair (A, B). If both a and b equal zero, then (A B) is singular
  210.      and S(I) = -1 is returned.
  211.  
  212.      An approximate error bound on the chordal distance between the i-th
  213.      computed generalized eigenvalue w and the corresponding exact eigenvalue
  214.      lambda is
  215.  
  216.           chord(w, lambda) <= EPS * norm(A, B) / S(I)
  217.  
  218.      where EPS is the machine precision.
  219.  
  220.      The reciprocal of the condition number DIF(i) of right eigenvector u and
  221.      left eigenvector v corresponding to the generalized eigenvalue w is
  222.      defined as follows:
  223.  
  224.      a) If the i-th eigenvalue w = (a,b) is real
  225.  
  226.         Suppose U and V are orthogonal transformations such that
  227.  
  228.                    U'*(A, B)*V  = (S, T) = ( a   *  ) ( b  *  )  1
  229.                                            ( 0  S22 ),( 0 T22 )  n-1
  230.                                              1  n-1     1 n-1
  231.  
  232.         Then the reciprocal condition number DIF(i) is
  233.  
  234.                    Difl((a, b), (S22, T22)) = sigma-min( Zl ),
  235.  
  236.         where sigma-min(Zl) denotes the smallest singular value of the
  237.         2(n-1)-by-2(n-1) matrix
  238.  
  239.             Zl = [ kron(a, In-1)  -kron(1, S22) ]
  240.                  [ kron(b, In-1)  -kron(1, T22) ] .
  241.  
  242.         Here In-1 is the identity matrix of size n-1. kron(X, Y) is the
  243.         Kronecker product between the matrices X and Y.
  244.  
  245.         Note that if the default method for computing DIF(i) is wanted
  246.         (see SLATDF), then the parameter DIFDRI (see below) should be
  247.         changed from 3 to 4 (routine SLATDF(IJOB = 2 will be used)).
  248.         See STGSYL for more details.
  249.  
  250.      b) If the i-th and (i+1)-th eigenvalues are complex conjugate pair,
  251.  
  252.         Suppose U and V are orthogonal transformations such that
  253.  
  254.                    U'*(A, B)*V = (S, T) = ( S11  *   ) ( T11  *  )  2
  255.                                           ( 0    S22 ),( 0    T22) n-2
  256.                                             2    n-2     2    n-2
  257.  
  258.  
  259.  
  260.  
  261.                                                                         PPPPaaaaggggeeee 4444
  262.  
  263.  
  264.  
  265.  
  266.  
  267.  
  268. SSSSTTTTGGGGSSSSNNNNAAAA((((3333SSSS))))                                                          SSSSTTTTGGGGSSSSNNNNAAAA((((3333SSSS))))
  269.  
  270.  
  271.  
  272.         and (S11, T11) corresponds to the complex conjugate eigenvalue
  273.         pair (w, conjg(w)). There exist unitary matrices U1 and V1 such
  274.         that
  275.  
  276.             U1'*S11*V1 = ( s11 s12 )   and U1'*T11*V1 = ( t11 t12 )
  277.                          (  0  s22 )                    (  0  t22 )
  278.  
  279.         where the generalized eigenvalues w = s11/t11 and
  280.         conjg(w) = s22/t22.
  281.  
  282.         Then the reciprocal condition number DIF(i) is bounded by
  283.  
  284.             min( d1, max( 1, |real(s11)/real(s22)| )*d2 )
  285.  
  286.         where, d1 = Difl((s11, t11), (s22, t22)) = sigma-min(Z1), where
  287.         Z1 is the complex 2-by-2 matrix
  288.  
  289.                  Z1 =  [ s11  -s22 ]
  290.                        [ t11  -t22 ],
  291.  
  292.         This is done by computing (using real arithmetic) the
  293.         roots of the characteristical polynomial det(Z1' * Z1 - lambda I),
  294.         where Z1' denotes the conjugate transpose of Z1 and det(X) denotes
  295.         the determinant of X.
  296.  
  297.         and d2 is an upper bound on Difl((S11, T11), (S22, T22)), i.e. an
  298.         upper bound on sigma-min(Z2), where Z2 is (2n-2)-by-(2n-2)
  299.  
  300.                  Z2 = [ kron(S11', In-2)  -kron(I2, S22) ]
  301.                       [ kron(T11', In-2)  -kron(I2, T22) ]
  302.  
  303.         Note that if the default method for computing DIF is wanted (see
  304.         SLATDF), then the parameter DIFDRI (see below) should be changed
  305.         from 3 to 4 (routine SLATDF(IJOB = 2 will be used)). See STGSYL
  306.         for more details.
  307.  
  308.      For each eigenvalue/vector specified by SELECT, DIF stores a Frobenius
  309.      norm-based estimate of Difl.
  310.  
  311.      An approximate error bound for the i-th computed eigenvector VL(i) or
  312.      VR(i) is given by
  313.  
  314.                 EPS * norm(A, B) / DIF(i).
  315.  
  316.      See ref. [2-3] for more details and further references.
  317.  
  318.      Based on contributions by
  319.         Bo Kagstrom and Peter Poromaa, Department of Computing Science,
  320.         Umea University, S-901 87 Umea, Sweden.
  321.  
  322.      References
  323.      ==========
  324.  
  325.  
  326.  
  327.                                                                         PPPPaaaaggggeeee 5555
  328.  
  329.  
  330.  
  331.  
  332.  
  333.  
  334. SSSSTTTTGGGGSSSSNNNNAAAA((((3333SSSS))))                                                          SSSSTTTTGGGGSSSSNNNNAAAA((((3333SSSS))))
  335.  
  336.  
  337.  
  338.      [1] B. Kagstrom; A Direct Method for Reordering Eigenvalues in the
  339.          Generalized Real Schur Form of a Regular Matrix Pair (A, B), in
  340.          M.S. Moonen et al (eds), Linear Algebra for Large Scale and
  341.          Real-Time Applications, Kluwer Academic Publ. 1993, pp 195-218.
  342.  
  343.      [2] B. Kagstrom and P. Poromaa; Computing Eigenspaces with Specified
  344.          Eigenvalues of a Regular Matrix Pair (A, B) and Condition
  345.          Estimation: Theory, Algorithms and Software,
  346.          Report UMINF - 94.04, Department of Computing Science, Umea
  347.          University, S-901 87 Umea, Sweden, 1994. Also as LAPACK Working
  348.          Note 87. To appear in Numerical Algorithms, 1996.
  349.  
  350.      [3] B. Kagstrom and P. Poromaa, LAPACK-Style Algorithms and Software
  351.          for Solving the Generalized Sylvester Equation and Estimating the
  352.          Separation between Regular Matrix Pairs, Report UMINF - 93.23,
  353.          Department of Computing Science, Umea University, S-901 87 Umea,
  354.          Sweden, December 1993, Revised April 1994, Also as LAPACK Working
  355.          Note 75.  To appear in ACM Trans. on Math. Software, Vol 22,
  356.          No 1, 1996.
  357.  
  358.  
  359. SSSSEEEEEEEE AAAALLLLSSSSOOOO
  360.      INTRO_LAPACK(3S), INTRO_SCSL(3S)
  361.  
  362.      This man page is available only online.
  363.  
  364.  
  365.  
  366.  
  367.  
  368.  
  369.  
  370.  
  371.  
  372.  
  373.  
  374.  
  375.  
  376.  
  377.  
  378.  
  379.  
  380.  
  381.  
  382.  
  383.  
  384.  
  385.  
  386.  
  387.  
  388.  
  389.  
  390.  
  391.  
  392.  
  393.                                                                         PPPPaaaaggggeeee 6666
  394.  
  395.  
  396.  
  397.